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("—^ I We report on our analysis of the Coulomb singularity problem in the frame of the coupled channel 

^vj , scattering theory including spin-orbit interaction. We assume that the coupling between the partial 

wave components involves orbital angular momenta such that Al — 0, ±2. In these conditions, the 
r^ ' two radial functions, components of a partial wave associated to two values of the angular momentum 

I, satisfy a system of two second-order ordinary differential equations. We examine the difficulties 

-^-s ' arising in the analysis of the behavior of the regular solutions near the origin because of this coupling. 

First, we demonstrate that for a singularity of the first kind in the potential, one of the solutions is 
not amenable to a power series expansion. The use of the Lippmann-Schwinger equations confirms 
this fact: a logarithmic divergence arises at the second iteration. To overcome this difficulty, we 

' fj ' introduce two auxilliary functions which, together with the two radial functions, satisfy a system of 

C/3 , four first-order differential equations. The reduction of the order of the differential system enables 

us to use a matrix-based approach, which generalizes the standard Frobenius method. We illustrate 
our analysis with numerical calculations of coupled scattering wave functions in a solid-state system. 
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I. INTRODUCTION 

The scattering of partial waves by a central perturbating potential is fully characterized by the scattering matrix. 
Beyond the Born approximation, this matrix, which reduces to a scalar in the absence of coupling, is simply expressed 
with the phase shift Si of the partial wave with orbital angular momentum I. When two waves with different orbital 
angular momenta are coupled by a spin-orbit interaction, the scattering matrix can be constructed with two amplitudes 



a^ and two phase shifts S^ as shown in Ref. H 

la i 



The generalized variable phase method@, yj is a powerful tool to obtain these four parameters: one defines two 
amplitude functions, a'^{r), and two phase functions, S^{r), which would characterize a scattered partial wave in a 
potential truncated at the radius r. These functions satisfy a system of four coupled non-linear first-order differential 
equations. The numerical solution of this system permits to obtain not only the four parameters in the asymptotic 
region (r greater than the range of the potential), but also the radial wave functions for each value of r. 

The generalized variable phase equations contain functions which are singular at the origin. As a consequence, 
application of standard numerical procedures right from the origin fail. In addition, the numerical problem is further 
complicated when the scattering potential is of the Coulombic type. Therefore, one is bound to analyze the behavior 
of the solutions of the differential equations near the origin. The computation of reliable solutions of the four coupled 
differential equations is otherwise impossible. 

The analysis presented in this article originates in our work on the scattering states of coupled semiconductor 
valence-band holes in a point defect potential [3|] ■ Our work now is devoted to the development of a general framework 
allowing the treatment of class of problems involving coupled Schrodinger equations of the form given by Eqs. ([2]) 
and ©. 

The article is organized as follows. In Section II, we introduce the definitions and notations we use throughout 
the paper. The mathematical framework is that of the standard time-independent nonrelativistic scattering theory 
presented in the books of, e.g., Newton[4], Taylorfsj, and Reed and Simon[(>]. In Section III, we see that for a singular 
potential only one of the regular solutions can be expanded as a power series; the search for an expansion of the 
other solution near the origin using the Lippmann-Schwinger equation leads to a logarithmically divergent integral. 
In Section IV, we develop a matrix-based method to treat this singularity problem. It yields an expansion of the 
radial wave functions near the origin, which may contain logarithmic terms. The presence of these terms explains 
the failure of the methods described in Section III. We provide numerical examples in Section V to illustrate our 
calculations. We choose the widely used screened Coulomb potential of the Yukawa type to model the scattering of 
holes by an ionized impurity potential Q, which corresponds to a realistic situation in doped semiconductors. Explicit 
expressions of functions and matrices omitted in the text, and a detailed discussion of a conditioning method for the 
Lippmann-Schwinger equation are given in a series of appendices. 

II. DEFINITIONS AND NOTATIONS 

We consider a Hamiltonian Ji that characterizes the scattering process of a particle with momentum Sfc, by a central 
potential V{r), where r is the distance from the origin of the potential. This Hamiltonian also contains a spin-orbit 
interaction term. We restrict our analysis to stationary scattering states of H with positive energy E. Introducing 
the variable x = kr, the partial waves characterized by the total angular momentum F = L + S, may be written 
asg: 

^=j:^^^^^^\L + e,S,F,F.), (1) 

where L and S are the orbital and spin angular quantum numbers, and Fz is the projection of F on the quantization 
axis. 

We assume that the radial Schrodinger equations satisfied by the functions ul±i{x) take the form: 



Bi 



d^ LiL-l) 



UL-l{x) +Bi2 



(P 2L + 1 d (L- 1)(L + 1) 
dx^ X dx x"^ 



UL + l{x) 



dx^ x^ 

= [v(x) - \\UL-1[X) (2) 



and 



dx2 


2i + l d L{L + 2)' 
X dx x'^ 


UL~l{x) +B22 



B21 



= [v{x) - l]uL+i{x) 



d^ (L + l)(L + 2) 
dx^ x"^ 



UL+l{x) 



(3) 



where v{x) is the reduced potential v{x) — V{x/k)/E. In the absence of the spin-orbit interaction, the matrix [B] 
whose elements are Bij {i,j = 1,2), reduces to the 2x2 identity matrix I2, and one recovers the uncoupled radial 
Schrodinger equations associated to the two values of the orbital angular momentum: L ± 1. The real matrix [B] is 
assumed to be symmetric and positive definite; it thus may be diagonalized as follows: 



[B] = [P][n]-^[P]-\ 



(4) 



where the matrix [T'] is a 2 x 2 passage matrix, whose elements are denoted Vij . The exact knowledge of the matrices 
[TZ] and [V] is not required to understand the analysis developped in the present article; we only assume that the 
diagonal elements of [TZ], Vp and ri„, are positive. 

III. BEHAVIOR NEAR THE ORIGIN: THE SERIES EXPANSION APPROACH 

The Schrodinger equations ^ and ^ can be recast as: 



[^] ^ + m [^] - [^] [^]) -i + i (I2 - [c]) [B] [c]]{llil^^^-Hx) - ^]{llil^0^ , (5) 



where [C] 



da;-^ 
L 



x dx x^ 
Hereafter, we assume that near the origin the behavior of the reduced potential v{x) 



-(L + 1 
is well represented by a convergent series expansion of the form: 



v{x) 



E 



VpxP, 



(6) 



where p is a positive integer and Vp are real coefficients. This amounts to restrict the nature of a possible singularity 
to a simp le p ole in a: = 0. Within this framework, unscreened and screened Coulomb potentials [14| . the square 
potential [15| as well as any potential with a linear behavior near the origin may thus be treated [i6| . The sum in the 
above definition of v{x) may, e.g., correspond to the contribution of screening in the Yukawa potential. 

A. Absence of coupling 

If we impose the equalities rp = r^ = 1, the matrix [B\ is equal to the identity matrix [see ^]. The differential 
system ([5]) thus reduces to: 



ii+i,.,Mq)M}(;;-g)=w.)-ii(;;-«) 



dx 



x^ 



(7) 



Given that the matrix product (I2 ~ [C\) [C\ is diagonal, the radial wave functions ul±i{x) are, in this case, solutions 
of two uncoupled Schrodinger equations: 



d^ , , /(/ + 1) , , . , M , X 

-—T ui(x) 5 — ui[x) + 1 - v{x)\ui(x) = 0, 

dx^ a;^ 



(8) 



with I = L±l. 



The free regular solution of (H)) is a Riccati-Bessel function^ u\ (x) = ji{x), which behaves like a;'+^/(2Z + 1)!! 
near the origin. The perturbed solution can thus be expanded as: 

oo 

uz(x)=a;'+i^a„a;" (9) 

near a; = 0. In ([9]), the coefficients a„ satisfy the recurrence relation: 

n 

(n + 1){21 + n + 2) a„+i = -a„_i + ^ V"p an-p, (10) 

with Oo arbitrary and oi = aoVo/{2l + 2). 

The series © can also be obtained iteratively from the Lippmann-Schwinger equation: 



' u\^\x) = u^^^x) + J^vix')gi{x,x')u\°\x')dx' 
^ u|")(a;) = u'i°\x) + J^ v{x')giix,x')u['''^\x')dx' 



(11) 



with 



gi{x,x') = ji{x')ni{x) ~ ji{x)hi{x'), (12) 

where fii is a Riccati- Neumann function[2]. With this method, the coefficients a„ are obtained by replacing the 

functions v{x'), gi{x,x') and u" {x') in (jlip by their expansions in ascending powers of x and x' [2,]. Note that at 
each step of the iteration the integrals converge at the origin. 

B. Presence of coupling 

1. Analysis of the series expansion 

If the coupling in the differential system ([5]) is restored, it seems a priori natural to adopt the same approach as 
above, i.e. to search for a power series of the wave functions ul±i{x). Let 



^ n— 1 ^ ^ 

In the above sum, there is no finite constant terms (ag, bg) to ensure the regularity of the wave functions: ul±i(x — 
0) = 0. From ([5]), we establish the following recurrence relation: 

int, + [£]) [B] (in + 1)1. - [£]) ( «-; ) = - ( J-J ) + g V^ ( J- ) (14) 

Note the similarity between ([T4|) and ([T0| . The matrix product (nl. + [C]) [B] is invertible except for n = L — 1 and 

n = L + 1. Using ([3T1) and the matrix [W], given in Appendix El we find that the regular free solutions m^^i of (O 
take the following form: 

ul^'_-^{x) = -a^rpViiJL-i{rpx) - a^rn^Vi2JL^i{rmx) (15) 

"^L+ii^) = o,^rpV2iJL+i{rpx) +aQr^V22jL+i{rmx) (16) 



where Gq and a^ are two arbitrary constants, which appear in the vector Cq (see Section IV. B). The two vahies 
of n for which the matrix product (nl2 + [C]) [B] is singular thus correspond to the orders of the Bessel functions 
appearing in the expression of the free solutions (fTS]) and (fT6| . 
By recurrence we obtain: 

• For n < L — 1: a„=0 and b„ — 



• For n = L — 1: 6l=0 and cl is an arbitrary constant. Note that by taking the lowest order in x of the right 
hand side of (IT5t we obtain: 



aL = -ia+rnr^+' + a,r^2ri+')/{2L 1)!! (17) 

• For n = L: 

• For n = L + 1: Vq 6i+i = and 

2SiiaL+2 + (2i + 3)Bi2bL+2 = [{Vi - 1) az. + Vq ul+i] /{2L + 1) (19) 

To satisfy the equation Vq &l+i = there are two possibilities: 

1. Vb = 0, in which case a^+i — and 6l+i — 0, and one of the constants 0^+2 or 6^+2 may be chosen arbitrarily 
to satisfy the following equation: 

2BiiaL+2 + {2L + 'i)Bi2bL+2 - (^^i - 1) aL/(2L + 1). 

With a choice of the constants oi, and 0^+2, two regular and independent solutions that can be expanded 
as a power series may thus be generated. 

2. Vo 7^ 0, which amounts to a Coulombic-type behavior in x = 0. In this case 6l+i = 0, which yields aL — Q and 
a^+i = 0. The equation 

26iiaL+2 + (2L + i)Bi2bL+2 = 

shows that only one arbitrary constant may be chosen and hence only one solution may be expanded as a 
power series. This solution may be obtained by cancelling the term in x^ in (|15p . which amounts to impose [cf. 



4^ii^P^' + ao"^i2'^m+' = 0. (20) 

At the lowest order this solution reads: 



--^(-) ^ - ^^ .^^ ^ r^.rl ( ,,,^-/^ ,, ) .-- + 0(.--) (21) 



UL+i{x) ) (2L+1)!! P 7512 '" P V -^ii/(2i + 3) 



2. Analysis of the Lippmann-Schwinger equation 

To understand why the sohition with ol ^ cannot be expanded as a power series, we study the Lippmann- 
Schwinger integral equations satisfied by the radial wave functions ul±i{x): 

UL-l{x) \ _ [ u^L-ii^) 1 



Ul+i(x) J \ u'-^l^ix) 



Ql-i.L-i(^,^') Ql-i.L+ii^.x')\ (uL-i{x') 



dx'v{x')Q{x-x')[ "^i-^-^r ' ro"''^+' ' ' ' (22) 

where the four functions Q^^i j^j^iix, x'), given in Appendix iBl define the Green's matrix of the unperturbed system, 
and Q{x) is the Heaviside step function. More precisely, we have to determine which type of function is generated by 
iteration from the free solution whose behavior is given by: 

At the first iteration, after cancellation of the terms in Vqx^^-^, Vix^ and V2X^~^^, the integral 



(0) Mi-'^M +^(-""^)- (23) 



j:enerates the term 



aL / w(a;')^°+i,L-i(a;,2;'))2;' dx' 
Jo 



^°"" '^'''^'' {rl - rl) .^+1 = 5,+ix^+i (24) 



AL{L + l) dot [V] 

in the function ul+i{x). Note that the calculations are rather tedious because of the contributions of the second order 
terms of the Bessel and Neumann functions ]l-i{x) and hL+i{x). 

At the second iteration, similar calculations show that the two integrals 

bL+i / v{x')gl±^L^-i^{x,x')x'^^^dx' 
Jo 

generate terms proportional to Vq in ul±i{x). Each of these terms contains one logarithmically divergent integral in 
zero: 



Notice that the above constant vector is colinear to the vector I , ^^ j in the solution (I?!]) , which can be expanded 

as a power series. Our analysis thus shows why only one solution can be expanded as a power series while logarithmic 
terms are expected in the other solution. 

It is now instructive to mention a conditioning method developped by Newton [J, [lO| . Considering the case of a 
tensorial scattering potential that couples two angular momenta of values I = L ± 1, Newton studied the behavior 
of the regular solutions near the origin and presented a method to obtain a well-conditioned Lippmann-Schwinger 
integral equation to avoid the divergence problem. A detailed description of Newton's method is given in Appendix 
[Cl Here, we can say that although Newton suggested that this method may be applied in a general case !4f], the 
modification of the Lippmann-Schwinger equation is impracticable in two particular situations: 

• when the difference between the values of the two coupled angular momenta is greater than 2; 

• when the coupling involves a nondiagonal Green's matrix. 

In the present work Al = 2 and the potential that couples the two components ul±i is scalar, but the Green's matrix 
is not diagonal and diverging terms appear only at the second iteration of the Lippmann-Schwinger equation. For 
this reason we have developed another method, presented in the next section. 



IV. BEHAVIOR NEAR THE ORIGIN: THE MATRIX APPROACH 
A. Reduction of the order of the differential system 

The coupled radial Schrodinger equations ([2]) and ((3)) may be written as a factor product: 

where A/ is a differential operator: A/ = d/da; + l/x. Having in mind that the matrix [B] can be diagonalized [Eq. 
(HI)], it is natural to introduce the two-component vector G{x): 

G{x) = [R]-^ [V]-' (^ ^0 ^ ^l^^ ) [V] u{x), (26) 

where the vector u{x) is defined by its components; Up{x) and Uni(a;), which are related to the radial wave functions 
ul±i{x) via: 

«p(^) \ [ppi (^ «L-i(^) y (27) 

Um{x) J \ UL + 1[X) J ' 

It follows that the second-order differential system ([25]) can be transformed ^ into a system of four first-order 
differential equations perturbed by a potential represented by a 4 x 4 matrix [t/(x)]: 

A(^|^^{[^(^)] + [T/(:,)]}(^|^, (28) 

The matrices [.4(x)] and [V^(a;)] are defined as follows: 

\A{x)\ = M (29) 

and 

({o\ {o\\ 

\V{x)\=v{x)\ , (30) 



l2 [O] 



where [O] is the 2x2 null matrix. 



B. Variable phase method 

The free solution of (1^5)1 . obtained for v(x) = for all x, is of the formQ 

^")-f'''^" '""^ imx)]Co, (31) 

where Co is a constant four-component vector. The 4x4 matrix [W(x)] contains the regular and irregular free 
solutions. It is given in the Appendix [Xj Since the differential system (P5|) is of order 1 , the Lagrange method of the 
variation of constants [9I| can be applied: 



[V]-' [O] . 

[W{x)] C{x), (32) 

[O] I2 



for v{x) 7^ 0. In this case, the 4-component vector C{x) is not constant. In the variable phase approach [2| it is 
searched as 0] 

(a+ (a:) cos (5+ (x) \ 
—a (x) sin (5 (.x) / 

where 5 {x) are the phase functions and a {x) are the amphtude functions. Note that for ease of notation, we omit 
the quantum numbers S, L, F, F^ in the equation above. 

Differentiation of the vector C{x) with respect to x yields a generalized form of the phase equations (GPE) [2[, 
which constitute a differential system, which we denote S (&=*=, 5^), satisfied by the functions i5^(x) and a^[x). Solving 
this system permits the computation of the wave function ul±i for all x; it also permits to find the four parameters o^ 
and 5^ involved in the scattering matrix. Indeed, these parameters are given by the asymptotic values: Mva^^^o or^{x) 
and Vmvx^oo 5^ {x) . The GPE must be solved with (5*(0) = as initial conditions to ensure the regularity of the 
radial wave functions ul±i{x) at the origin. 

The differential system S(a^, J ) contains the irregular functions hi. In some cases, the reduced potential, v{x)^ 
may diverge too when x — > 0. The computation of the solutions of the GPE thus requires an in-depth analysis of the 
behavior of the functions a^{x) and 5^{x) near the origin, which can be deduced from those of u and G [see Eq. (|32p ] 
using a generalized Frobenius method. 

C. Generalized Frobenius method 

Near the origin, the matrix [yl(a;)] + [^(a;)] in (pS)) can take the form of an expansion: 

00 
[A{x)] + [V{x)]^Y.'^A,]xP-\ (34) 

where the matrices \Ap\ are given in the Appendix[Dl Since the matrix [^0] is never zero, &x~^ term, called singularity 
of the first kind, appears in the expansion of the matrix [.A(x)] + [y (a)] near the origin. 
The general solution of (P5|) near the origin can be searched for as [l7| : 

U \ „ v-^ v-^ T» X-' 



G 



Y^Y^C,., -[Hx)]\ (35) 



where i and j are integers that satisfy: < z < 3 for a 4 x 4 system and j > 0, t] is an eigenvalue of the constant 
matrix [Aq]. The series expansion in psp may contain terms in ln(a;), [ln(x)] and [ln(a;)] . Using (|28p. it can be 
shown that the four-component constant vectors C^.j satisfy the following recurrence relation: 



{{J + V)U - [Aq]}C^,j + C,+ij = Y ["^"] <^^J-"' (36) 

ri=l 

with the condition that for all j, C4J = 0. One method to obtain systematically the relevant powers of x and 
ln(x) is given in [17J. In the present work, the problem is reduced to the diagonalization and inversion of 4x4 
matrices. One can check that the four eigenvalues of the matrix [Aq] are those of the matrices [£] and — [£], i.e. 
{—L — 1; —L; L; L + 1}. These eigenvalues correspond to the power of a; of the leading term of the expansion of the 
special functions zi (x) — hi (x) near x = or zi (x) = ji (x) [18| , as shown in Table HI 



TABLE I. Values taken by the power 77 of the leading term of the expansions of the special (Riccati-Bessel or Riccati-Neumann) 
functions zi(x) near a; = 0. 



zi{x) 



-L-1 flL+lix) 

-L nL{x) 

L JL-i{x) 

L + 1 Jl{x) 



D. Properties of the vectors Cij{a,b) associated to the regular solutions 

Solving the two coupled second order differential equations, psp . satisfied by the functions ui±i{x), yields four 
constants of integration that correspond to the four independent solutions, of which two are regular and two are 
singular. The two smallest eigenvalues correspond to two singular solutions which are disregarded by determining the 
vectors Cij so that: Ci_j = for < i < 3 and < j < 2L, which amounts to impose that the related two constants 
of integration are zero. The choice of the values of the other two constants corresponding to the regular solutions 
remains free. The consequences of this choice are now discussed. 

1. Eigenvalue rj — L 

This eigenvalue corresponds to a behavior in x^, which is acceptable. The vector Cij is now associated to the term 
x^+^ [ln(a;)]\ Let us consider the case j — for which the smallest power of x in the expansion of (|35l) coincides 
with the eigenvalue rj = L. Given that the coefficients Ctj of the irregular solutions are equal to zero, the recurrence 
equation, (1551) . reduces to: 



for i = 0, 1, 2, 3. Since C4 = 0, the above equality yields: 



(37) 



{il4 - [^0]} Co,o = 0, (38) 

which implies that either Co,o = or Co,o is an eigenvector of the matrix [^0] associated to the eigenvalue r/ — L: 



^2L[Vr'(l^^ 



Co,o(a) = a 



Vo\V] 







(39) 



/ 



where a is an arbitrary constant. In both cases, we also find that Ci.o = C2,o = C3.0 = 0, which is equivalent to the 
absence of x^ [ln(x)]' terms and corresponds to the presence of the Riccati-Bessel function of order i — 1 in ul^i{x). 



2. Eigenvalue r] — L + 1 



For J = 1, the recurrence relation reads: 



{{L + l)l4 - [A]}C^,i +C,+i,i = [A]C,,o, 
from which we deduce the equation satisfied by the vector Co,i: 



(40) 



{{L + l)l4 - [A]}V"o,i = {{L + l)l4 - [A^]f [A]Co,o. 



10 



(41) 



Since the matrix {L + l)l4 — [^o] is not invertible, the above inhomogeneous matrix equation has a solution 
only on the condition that the vector on the right hand side of (HIT) is orthogonal to the kernel of the operator 
{(1/ + 1)14 — [^0]} [13 • One can check that this condition is fulfilled when Co,i takes the form: 



( Vo[v]-' [v]-^ [v][nf [T]-' [l^ ^ 



Co,i(a, &) = a 



\ 



^o[V] 



( 



\ [^r^ \ 



\ 



(42) 



where the matrix \D] and the constant ^0 are given by: 



[V] = 



1 







2L + 2 



(43) 



and 



Co- 



1 



2i + l 



2L{Vi - 1) + Vc 



^ViiV22rl-Vi2V2irl, 



det [V] 



(44) 



The constant a is the same that appears in the definition of Co,o(a)j (|39p . and h is the other arbitrary constant that 
has to be introduced to obtain the general solution. These two constants are indeed necessary to generate the two 
independent solutions that are regular at the origin. Once the vectors C^.o and Co,i are known, it is possible to 
calculate the vectors Ci^i for 1 < i < 3. Note that the second vector on the right hand side of (UU is an eigenvector 
of the matrix [^0] associated to the eigenvalue 77 = L + 1; as such it does not bring any contribution to the vectors 
Ci,i for 1 < z < 3. Indeed we find: C2,i = Cs,! = since Ci,i is also an eigenvector of the matrix [^0] associated to 
1] = L + 1, and 



/ 



Ci,i(a) = qVq 



a [v]-' 



\ 



1;/ 



(45) 



where ^1 is given by: 



1 ^21^22 {rl - rl) 



2L + 2 



det [V] 



(46) 



Note that the last two components of the vector Ci,i give rise to terms in x^+^ [ln(a;)], which are of the second order 
in the potential strength. 

3. Case j >2 

Since from J = 2 the matrices (L + j)l4 — [.4o] are invertible, the vectors Cij can be determined without the need 
to introduce new constants. For j — 2, we obtain 



Ci,2(a) = {{L + 2)l4 - [^0]} [A] Ci,i(a) 

^2,2 = Cz,2 = 

Co,2(a, h) = {{L + 2)l4 - {Aq\Y^ { [^1] Co4(a, h) + [^2] C^M - C^M\ 



(47) 
(48) 
(49) 



11 



It is not useful to give the explicit form of the vector Ci_2, but we stress that a lengthy calculation yields four 
non-zero components. Terms in x^~^'^\n{x) thus appear in the wave functions Up{x) and u^{x), but they do not 
compromise their regularity at the origin. 

To check if terms in [ln(a;)] and [ln(a;)] do occur in the series ([35]), we consider the general equations for j > 2: 



{ (L + j)l4 - [Ao]}Co,j + Cij = J2 [-^»] ^0,j-n (50) 

n=l 

{{L + j)U - [Ao]} Ci,j + C2j = J2 i-^n] Cl j-n (51) 

n=l 

{{L + .7)14 - [Ao]}C2,j + Csj = J2 [-^«] C"2j-n (52) 

n=l 

{{L + j)U - [Ao]}C3.j - Y. [-^"1 ^"sj-" (53) 

n=l 

Recalling that Ci^ — €2,0 = Ca^o — C24 = Cs^i = 0, we find by reccurence C2J = Csj- — 0. The series (|35l) thus does 
not contain terms neither in [ln(x)] nor in [ln(x)] . This result allows us to express the vector Cij as a function of 



the vectors Ci,fc, k < j — 1: 



Cij(a) = {{L + j)U - [A^]} ' Ij2 [A^]C;,,_„(a) j , 



_ (54) 

\n=l / 

which we use to obtain the vector Cq.j as: 

4,(a,&) = {(L+j)l4-[A]}'' (^[X]Co.j-n(a,&)-Ci,,(a)j . (55) 

We observe that the vectors Cij associated to the terms in x^^^ ln(a:), are proportional to the factor aVQ{r^ — r^) 

appearing in the expression of the vector Ci^i, (|45l). The ln(a;) terms are thus only generated if the potentials contain 
a Coulombic component (Vb ^ 0) and if the two differential equations satisfied by the wave functions Up{x) and Um{x) 
are coupled (r^ 7^ ''p)- Since this factor is of second order in the potential strength, we find that it is consistent with 
the fact that the divergence problems appear only from the second iteration when searching for the solutions using 
the Green's matrix Q. 

E. Logarithm-free solution 

The choice a = leads to Co,o = and Ci,j — for all j. As a consequence, we see that there exists a solution 
whose expansion does not contain any logarithmic term. This particular expansion is called a Frobenius series |13l , |2l| 
and here it takes the following form: 

^) =bx^+^Y.'^-'o,,+i{0,l)x^, (56) 

where the vectors Co.j, j > 2, now satisfy the following recurrence relation: 

Co J - {{L+j)U - [Ao]}-' IYI [■^n]Co,j-n] ■ (57) 
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The terms of the lowest order of the functions Up{x) and u^{x) are in x^"*"^ since the two first components of the 
vector Co,i are zero [see dH])]; more precisely, we obtain: 



1 
2i + 3 

at the lowest order of the expansion. 



h\^ \[V][nf[V]-'[ U^+2 (58) 



F. Structure of the general solution 

From the definitions of the vectors Co.o(a) and Co,i(a, &), and from the linearity of the recurrence relations satisfied 
by the vectors Co,j and Ci.j, the following properties can be deduced: 

Co,, (a, &) = a Co,,(l, 0) + h Co.,(0, 1) (59) 

C;,,(a)=aCi,,(l) (60) 

The vectors Co,j(0, 1) and Co,j(l, 0) are determined uniquely, recursively, from the series expansion of [.4(x)] in a; = 0. 

Moreover the vectors Cij(l) and Co.j(0, 1) satisfy the same recurrence relation. Since the first non-zero vectors satisfy 
the following equation: 

^^■^(^^ - "2LT^ ^^\ ^°'^(°' ^^ ^^^^ 

we obtain 

2 -P„,'P„„/'^2 _ 2\ 



V ^2lT'22 (r 



T^ 



so that the general regular solution generated by two arbitrary constants (a, b) can be written [20[: 

^/(a,6) V<^/(1,0) V<^/(0,1) 

where 



u 

(0,1) j>0 



g; =5]4,(0,1)x^+^ (64) 



and 



= 5]Co,,(l,0).T^+^ 
(1,0) j>0 



aVp^ ^21^22 (rg, - rg) /^ 
2L + 2 det [75] V G 



(0,1). 



ln(a;). (65) 



The vectors =; j and ^ 1 are thus two independent solutions, of which the first is the logarithm-free 

V'-^/(o,i) V'^/(i,o) 

solution. 

If a ^ 0, the expansion consists of a power series plus an additional term originating from the product of Vq \n{x) 
and the previous solution obtained with (a = 0, 6 = 1). This ln(a;) term, studied in the theory of Fuschian differential 
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equations [i2,ll3 appears only at the second order for Coulombic potentials, i.e. v{x) w Vox~^ when a- ^- 0. This is 
consistent with the results obtained from the Lippmann-Schwinger equation. 
The relationships between the constants (a, 6) and (aq ,aQ") are as follows: 



^-I^o\ 


[v] [7^]^+' 


1 ) 


(2L + 1)!! 



(66) 

'0 / 

They are obtained from the comparison of the lowest order terms of the series expansion of uu.h) in (1631) . and u^j_i{x) 
in P^ and ([TB| . One can check that for the logarithm- free solution the equality (pn|) is satisfied. 

V. NUMERICAL EXAMPLE 

From a computational viewpoint, the power series appearing in the expansion of (^ ) must be truncated, which 

limits its use to a finite interval [0;Xconv]; moreover Xconv is in general smaller than the distance where the functions 
a^{x) and S^{x) cease to vary and which is directly linked to the range of the potential v{x). To obtain a^(oo) and 
(5 (oo), the differential system S(a , S ) must be numerically solved from a value of a;, a; = x^ < Xconv The values 
of a^(xm) and S {xm) at the matching point x = Xm are found from Eqs. (l32t and (|33|) . 

To illustrate our work we choose the case of the spherical Hamiltonian of semiconductor valence-band holes in a 
perturbating spherical potential V{r), which reads [8|: 

n = ^p'-^ P(2) • 5(2) + V{r), (67) 

2mo 6mo 



where mo is the free electron mass, and 71 and fi are two material-dependent parameters. The first term of "H in (|67|) 
is the kinetic energy. The second term represents the spin-orbit interaction. The operator P^^> ■ S''-^' is the scalar 
product of two irreducible tensors of rank 2 constructed from the components of the hole momentum, p, and the spin 
angular momentum, S [8|. 

The quantum numbers of the scattered partial wave are chosen as i^ = 2 and L — 2. With this choice, the matrix 
[B] is readily obtained [8,]: 



The matrices [V] and [7^] easily follow: 



l^l'-Ks^f'f) (^«' 



pi-;^(":l:l)' <«" 



and 



The interaction between the ionized defect and the holes is represented by a Yukawa potential: 

, , 2Z 6"^/'=^= , , 

«(^) = I ' (71) 

/cob X 

where Z is the charge number, ob is the Bohr radius[3|, and Ag the screening length parameter. This corresponds to 
the Thomas-Fermi approximation for the screening. We took the following values: k = 2/aB, Z = —1, Xg = 20aB, 
and /i = 0.481 for the effective mass parameter of Si \2^. The differential system T.(a^, 6^) for this case is similar to 
that given in reference [3J. 
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FIG. 1. Scattering amplitudes a (x) and a^{x) as functions of the scaled distance x = kr. 
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FIG. 2. Scattering phase shifts 5 (x) and 5'^{x) as functions of the scaled distance x = kr 



A. Logarithm-free solution 



As discussed in Sec. IV. E, the logarithm-free solution is obtained setting a = 0. The amplitude, phase and radial 



and 



10Pi2/r^, for 



wave functions are shown on Figs. 1, 2, and 3. This can be done by taking Cg = — lOPn/rj^ 
L = 2 [see (gni)]- 

On figure [H we see that the constants Oq and Og indeed correspond to the initial values of the amplitudes. The 
dashed-point curves on the three figures are obtained from the expansion series, ((55)) . limited to iV = 30 terms. 
Beyond the matching point Xm (here ccm = 3), the solid lines correspond to the solutions obtained by solving the 
differential system satisfied by a (x) and 6 [x) with the Runge-Kutta method ^]. Some experimentation is needed 
to find the minimum value of N for which the curves generated from the expansion remain unchanged in the interval 
[0;a;m], and match with those obtained by the Runge-Kutta method. Increasing the value of N may only allow to 
obtain a matching over a wider range; moreover, this presents little interest since this requires more computing time. 
At low energy (fcae <C 1), the oscillatory behavior of the solutions near the origin is enchanced [3| and the value of 
Xm must be reduced accordingly. 



B. Solutions with logarithm 



The solutions with logarithm were generated taking Og = a^ = 1. The curves shown on Figs. 4, 5 and 6 represent 
the amplitude, phase and wave functions. We obtain a good agreement between the series expansion (dashed lines) 
and the Runge-Kutta method (solid lines) for Xm = 3 and N — 30. In this case, a^ and aj do not correspond to the 
initial values of the amplitude functions Og (x), see figured 

Note that the range over which we performed the numerical calculations needs to be increased to obtain Ralph's 
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FIG. 3. Radial wave functions Up{x) and U3{x) as functions of the scaled distance x = kr. 
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FIG. 4. Scattering amplitudes a (x) and a'^{x) as functions of the scaled distance x — kr. 



parameters [ij, a^(oo) and 5^{oo). 



VI. CONCLUSION 

To compute the wave functions of the scattering states in a Coulombic potential, one must treat the related 
divergence problems with special care in the presence of coupling. We showed that only one of the two regular 
solutions can be expanded as a power series contrary to cases treated in the literature [Ifl, lll|- We presented a 
matrix-based method that makes it possible to seek the expansion of the two solutions in a systematic way. Using 
this method, we showed that for one solution the expansion near the origin contains logarithmic terms. These terms 
appear only at the second iteration of the Born series. This behavior, which is related to the particular form the 
spin-orbit interaction, explains why a well-conditioned Lipmmann-Schwinger integral equation is extremely difficult 
if not impossible to obtain. To our knowledge the solutions of the spin-orbit-coupled Schrodinger equations had not 
been studied before by direct numerical computation. 
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Appendix A: Matrix [W] 

For convenience, the 4x4 matrix [W] can be written as: 



[W] = 



[W^ii] [W12] 



[W21] [W22] 

where the four 2x2 blocks [Wn], [W^i2], [VF21], and [VF22] are given by: 



[W^ii] = r-p 



\Wl2] - rn 



-Vii]L-i{r^x) -ViifiL-iirpx) 
'P2iJL+i{rpx) 'P2inL+i{rpx) 



-Vi2JL-i{rmx) -VufiL-iir.j^x) 

'P22]L+l{rmX) 'P22nL+l{rmX) 



[W2 



hirpx) fiLirpx) 











(Al) 



(A2) 



(A3) 
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[W22] = ( ) (A4) 



Appendix B: Green's matrix elements 

The expressions of the four matrix elements of the Green's matrix [(y°] defined in (|22|) are given befow: 






nO f '\ ^p' 11' 22 , /N 'm/-^12/-^21 / /n ,r3T\ 



,0 / i\ rp'Pi2'P2i , ,s , raJ>ii'P22 



^0 / /N 'p/-^12/-21 , ,^ , 'm/-^ll/-^22 , ,-. .„„^ 



Gl-ul+i{x, x') ^ ^^^ ^ V^gL-ix+ArpX, r^x') - rm5L-i,L+i(r-ma;, r^a;')] (B3) 

^L+i,L-i(2;, x') = ^^^ ^ h?'p5L+i,L-i(''pa;, rpx') + r„igL+i_L_i(rijiX, r^x')\ , (B4) 



where 



gi,i'{x,x') = -j/(a;)n;'(a;') + n;(a;)jp(a;'), 
which is a form more general than (|12p . 



Appendix C: Newton's conditioning method 

By putting together the two components ^]°_;^(a;) and *}^^j(x) of the two independent regular solutions (a — a, b) in 
the same matrix [^(cc)], the Lippmann-Schwinger integral equation corresponding to the tensorial scattering potential, 
takes the following form: 

f^i'Ux) ^tM\ _ f JlMx) 



\^i%{x) ^^^lix) ) V h+iix) 

Note that the Green's matrix, whose elements are the functions gi{x,x') defined in (IT^ . is diagonal in this case. The 
first term on the right hand side of (|C1[) is the 2x2 matrix [^'^"^(x)] , solution of the equation at the order 0. The 

behavior of the functions ^}J^j^(x) near the origin a; = is analyzed by iteration. For the solution (b) the integrals 

converge at all orders and the functions 4'}^^-^(x) can be expanded as power series. For the solution (a) the first 
iteration of the Born series leads to the evaluation of the following integral: 

v{x')gL+iix,x')JL-i{x') dx', 







which usually diverges because of the presence of the product nL+i{x')JL-i{x') that behaves like — (2L + l)/x' near 
x' = 0. We therefore see that problems arise when the index I of the Green's function gi is greater than that of the 
Bessel function J;/ , which represents the free solution. 
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These difficulties are circumvented by adding a term to the Lippmann-Schwinger equation, (IC1[) . which modifies 
the inhomogeneity in a convenient wav|4l.[l0|: [\I>°(a;)] becomes [^°(x)] (I2 + [M]), where [AI] is a constant matrix 
appropriately chosen so that it cancels the problematic term that reads 

^(2L + l)ii+i(a;)£t;(a;')(i/^, g) dx' 
The appropriate matrix [M] for the case considered, may be written as follows: 

[M] = ~{2L + 1) / - [yOFF(^,)] d^,^ (C2) 

Jo 2; 

where [T^'^^^(a;)] = 1 / \ n ) ' ^^'^ a^o > may be arbitrarily chosen. 

The integral that appears in the definition of [M] in (IC3p may be decomposed into two parts as: /„ "(. . .)dx' = 
L (. . .)da::' + J °(. . .)da;', to obtain the modified Lippmann-Schwinger equation: 

[nx)]^[^^'H-)]{t,-{2L + l)£^ [V^^^ix')]^ (C3) 

+ J^^ (^[g'{x,x')] [V{x')] [vl/(x')] - ^^ [l^°^^(x')]) dx', 

This equation makes it possible to obtain equivalents to the solutions (a) and (b) in the vicinity of x = 0. Note that 
it is no longer necessary to assume that the integral that defines the matrix [M] converges. 

It also is instructive to see how the solutions (a) and (b) behave when the tensorial potential of (jCip is constant 
|11| . Taking xq = 1, we obtain at the lowest order, a modified solution of the form: 

'¥°Hx)\il2 + [M])=( f^,^,.^''-}^''\ ^, . . ^ ^ f ^) (C4) 

L ^ 'J ^ \{2L + l) V jL+i[x)\n{x) jL+i[x) J 

where v is the off-diagonal term. The second component of the solution (a) contains an amplitude function oj^Ax) 
that may diverge logarithmically. This behavior, predicted in the theory of Fuschian differential equation |12l. Il3|. 
does not compromise the regularity of [$] in x = 0, since Jl+i{x) w x^~^^. 

Appendix D: Matrices [^0]) [-4i] and [An] 



[V] ^[C][P] [O] 

{Ao]={ ^ I, (Dl) 

[o] [nf 

[Ai] =1 I , (D2) 

(Fl-l)l2 [O] 

where [O] is a 2x2 matrix whose elements are all zero; and for n > 2 

f[0] [0]\ 
[An] = K . (D3) 

V I2 [o] 
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